Measuring schroeder periods
Fine scale acoustic perception in zebra finches
Source code and data found at https://github.com/maRce10/acoustic-fine-features-zebra-finch
1 Purpose
- Detect start and end of periods in a schroeder using amplitude cross-correlation
- Compare methods for measuring fine scale structural variation and periodicity in schroeders
2 Report overview
3 Synthetizing schroeders
3.1 Function to make schroeders
Code
qwindpc <- function(rftime, srate, sig) {
ramppts <- round((rftime/1000) * srate)
hold1 <- seq(0, 1, length.out = ramppts + 1)
onramp <- sin(0.5 * pi * hold1)^2
offramp <- cos(0.5 * pi * hold1)^2
steady <- rep(1, length(sig) - (2 * ramppts) - 2)
wind <- c(onramp, steady, offramp)
if (length(wind) > length(sig)) {
print("Window not applied, block too short for window")
} else {
sig <- wind * sig
}
return(sig)
}
make_schroeder <- function(samp.rate = 44100, f0, dur = 1000, save.wave = FALSE,
plot = TRUE, n.components, scalar = 1, color = wave_col, path = ".",
file.name = NULL, random.start = FALSE, mfrow = c(2, 1), par = TRUE,
cex.oscillo = 1, ...) {
# duplicate duration if starting randomly
dur.mult <- if (random.start)
2 else 1
# set time instances at which samples will be 'taken'
t <- seq(1/samp.rate, 1, by = 1/samp.rate)
sumwavepos <- rep(0, samp.rate)
durpts <- round((dur * dur.mult/1000) * samp.rate)
phase <- rep(0, n.components)
amplin <- rep(1, n.components)
f <- rep(0, n.components)
Ns <- 1:n.components
for (n in Ns) {
compnum <- n
phase[compnum] <- scalar * pi * compnum * (compnum - 1)/n.components
f[n] <- f0 + ((compnum - 1) * f0)
}
wave <- matrix(0, nrow = n.components, ncol = length(t))
for (i in Ns) {
wave[i, ] <- amplin[i] * cos((2 * pi * f[i] * t) + phase[i])
sumwavepos <- sumwavepos + wave[i, ]
}
# make sumwapos longer so it fit the target duration
org.sumwavepos <- sumwavepos
while (length(sumwavepos) < durpts) sumwavepos <- c(sumwavepos,
org.sumwavepos)
posschr <- sumwavepos[1:durpts]/max(abs(sumwavepos[1:durpts]))
# plot(t[1:length(posschr)], posschr, type = 'l')
spectr <- fft(sumwavepos)
power <- abs(spectr)^2
sample <- length(t)/max(t)
freq <- (1:length(t)) * sample/length(t)
powerdb <- 10 * log10(power[1:(length(t)/2)])
if (random.start) {
strt <- sample(1:(length(posschr)/2), 1)
posschr <- posschr[strt:(strt + (length(posschr)/2) - 1)]
}
posschr <- qwindpc(10, samp.rate, posschr)
posschr <- posschr/max(abs(posschr))
wave_obj <- tuneR::normalize(tuneR::Wave(posschr, samp.rate = samp.rate,
bit = 16), unit = "16")
if (plot) {
opar <- par()
opar$cin <- opar$cra <- opar$sci <- opar$cxy <- opar$din <- opar$csi <- opar$page <- NULL
# on.exit(par(opar))
if (par)
par(mfrow = mfrow, mar = c(5, 4, 0, 0) + 0.1)
plot(freq[1:(length(t)/2)], powerdb, type = "l", xlim = c(0,
5200), xlab = "Frequency (Hz)", ylab = "Power(dB)", col = color,
...)
seewave::oscillo(wave = wave_obj, colwave = color, cexlab = cex.oscillo,
cexaxis = cex.oscillo)
}
# Save audio to a file
if (save.wave) {
if (is.null(file.name)) {
file.name <- paste0("f0-", f0, "_ncomp-", n.components,
"_", if (scalar == 1)
"pos" else "neg", ".wav")
}
writeWave(object = wave_obj, filename = file.path(path, file.name),
extensible = FALSE)
} else {
return(wave_obj)
}
}3.1.1 Schroeder at 700 Hz with 7 components (harmonics)
Positive
Code
Negative
3.1.2 Schroeder with 10 components
Positive
Code
Negative
Code
Randomly starting at any point
Code
Make example video:
Code
ani.options(ani.res = 200, ani.width = 1700, ani.height = 1000, ani.type = "png")
## make sure ImageMagick has been installed in your system gifs
## <- function(lng_dat, type = 'gif'){
# anni_fun <- if(type == 'gif') saveGIF else saveVideo
saveGIF(expr = {
for (i in 1:13) {
par(mfcol = c(2, 2), mar = c(5, 5, 3, 2) + 0.1)
ms <- make_schroeder(n.components = i, samp.rate = 44100,
f0 = 400, dur = 25, save.wave = FALSE, plot = TRUE, scalar = 1,
color = viridis(10)[3], random.start = FALSE, par = FALSE,
main = paste0("Positive, ", i, " components"), cex.main = 2,
cex.lab = 2, cex.oscillo = 2)
ms <- make_schroeder(n.components = i, samp.rate = 44100,
f0 = 400, dur = 25, save.wave = FALSE, plot = TRUE, scalar = -1,
color = viridis(10)[8], random.start = FALSE, par = FALSE,
main = paste0("Negative, ", i, " components"), cex.main = 2,
cex.lab = 2, cex.oscillo = 2)
}
}, video.name = file.path("./output/movies", "schroeders_by_sign_and_harmonics2"),
interval = 1)
# }, movie.name = paste0(z, '.', type), interval = 0.05) }
3.2 Make schroeders
Full factorial design of the following features:
- Components: from 1 to 25 (21 values)
- Fundamental frequency: from 200 to 800 (7 values)
- Sign: 1 and -1
- 294 schroeders
Create schroeders as individual sound files:
- Fixed start
Code
shr_tab <- expand.grid(type = c(-1, 1), components = seq(5, 25, 1),
f0 = seq(200, 800, 100))
nrow(shr_tab)
shr_tab$name <- paste0(f_pad_zero(x = 1:nrow(shr_tab), width = 3),
"_f0-", shr_tab$f0, "_ncomp-", shr_tab$components, "_", ifelse(shr_tab$type ==
1, "pos", "neg"), ".wav")
out <- warbleR:::pblapply_wrblr_int(seq_len(nrow(shr_tab)), function(x) {
make_schroeder(n.components = shr_tab$components[x], samp.rate = 44100,
f0 = shr_tab$f0[x], dur = 1000, save.wave = TRUE, plot = FALSE,
scalar = shr_tab$type[x], color = "red", path = "./data/processed/schroeders",
file.name = shr_tab$name[x])
})- Random start
Code
shr_tab <- expand.grid(type = c(-1, 1), components = seq(5, 25, 1),
f0 = seq(200, 800, 100))
nrow(shr_tab)
shr_tab$name <- paste0(f_pad_zero(x = 1:nrow(shr_tab), width = 3),
"_f0-", shr_tab$f0, "_ncomp-", shr_tab$components, "_", ifelse(shr_tab$type ==
1, "pos", "neg"), ".wav")
out <- warbleR:::pblapply_wrblr_int(seq_len(nrow(shr_tab)), function(x) {
make_schroeder(n.components = shr_tab$components[x], samp.rate = 44100,
f0 = shr_tab$f0[x], dur = 1000, save.wave = TRUE, plot = FALSE,
scalar = shr_tab$type[x], color = "red", path = "./data/processed/random_start_schroeders",
file.name = shr_tab$name[x], random.start = TRUE)
})Put Schroeders into a single sound file
Code
# fix start
est_schr <- warbleR::selection_table(whole.recs = TRUE, path = "./data/processed/schroeders",
extended = TRUE, confirm.extended = FALSE)
est_schr$bottom.freq <- (as.numeric(gsub("f0-", "", (sapply(strsplit(est_schr$sound.files,
split = "_"), "[[", 2)))) - 50)/1000
est_schr$top.freq <- freq_range(X = est_schr, threshold = 2, fsmooth = 0.8,
parallel = 13)$top.freq
est_schr$f0 <- gsub("f0-", "", sapply(strsplit(est_schr$sound.files,
split = "_"), "[[", 2))
est_schr$label <- gsub("os.wav_1|eg.wav_1", "", sapply(strsplit(est_schr$sound.files,
split = "ncomp-"), "[[", 2))
saveRDS(est_schr, "./data/processed/extended_selection_table_schroeders.RDS")
master <- master_sound_file(X = est_schr, file.name = "schroeder_master",
dest.path = "./data/processed/", gap.duration = 0.8, cex = 14)
write.csv(master, "./data/processed/master_annotations_schroeders.csv",
row.names = FALSE)
master$f0 <- c("start_marker", paste("f0", est_schr$f0), "")
master$label <- c("marker", est_schr$label, "marker")
Rraven::exp_raven(master, path = "./data/processed/", sound.file.path = "./data/processed/")
warbleR::full_spectrograms(master, path = "./data/processed/", sxrow = duration(readWave("./data/processed/schroeder_master.wav"))/12,
rows = 12, horizontal = TRUE, dest.path = "./data/processed/spectrograms",
collevels = seq(-150, 0, 5), labels = "label", song = "f0", wl = 2000)
# random start
est_schr <- warbleR::selection_table(whole.recs = TRUE, path = "./data/processed/random_start_schroeders",
extended = TRUE, confirm.extended = FALSE)
est_schr$bottom.freq <- (as.numeric(gsub("f0-", "", (sapply(strsplit(est_schr$sound.files,
split = "_"), "[[", 2)))) - 50)/1000
est_schr$top.freq <- freq_range(X = est_schr, threshold = 2, fsmooth = 0.8,
parallel = 13)$top.freq
est_schr$f0 <- gsub("f0-", "", sapply(strsplit(est_schr$sound.files,
split = "_"), "[[", 2))
est_schr$label <- gsub("os.wav_1|eg.wav_1", "", sapply(strsplit(est_schr$sound.files,
split = "ncomp-"), "[[", 2))
saveRDS(est_schr, "./data/processed/extended_selection_table_random_start_schroeders.RDS")
master <- master_sound_file(X = est_schr, file.name = "random_start_schroeder_master",
dest.path = "./data/processed/", gap.duration = 0.8, cex = 14)
write.csv(master, "./data/processed/master_annotations_random_start_schroeders.csv",
row.names = FALSE)
master$f0 <- c("start_marker", paste("f0", est_schr$f0), "")
master$label <- c("marker", est_schr$label, "marker")
Rraven::exp_raven(master, file.name = "random_start_schroeders", path = "./data/processed/",
sound.file.path = "./data/processed/")
4 Detecting periodicity
- Create a function that returns each detected segment in a list
- Two methods:
- Empirical Mode Decomposition (EMD)
- Time Autocorrelation (ac)
- The function can plot the mean period +/- standard deviation
Code
est_schr <- readRDS("./data/processed/extended_selection_table_schroeders.RDS")
mean_segment <- function(wave, cores = 1, plot = TRUE, pb = TRUE,
thinning = 1, col = wave_col, mean = TRUE, type = "ac", npeak = 20) {
# thin
if (thinning < 1) {
if (length(wave@left) * thinning < 10) {
stop2("thinning is too high, no enough samples left for at least 1 sound file")
}
# reduce size of envelope
wavefrm <- stats::approx(x = seq(0, duration(wave), length.out = length(wave@left)),
y = wave@left, n = round(length(wave@left) * thinning),
method = "linear")$y
} else {
wavefrm <- wave@left
}
# get empirical mode decomposition
if (type == "EMD") {
emds <- EMD::emd(wavefrm, seq_len(length(wavefrm)), boundary = "wave")
perd <- emds$imf[, 4]/max(emds$imf[, 4])
# plot(x = seq_len(length(wavefrm)), y = perd, type =
# 'l') lines(y = wavefrm / max(wavefrm), x =
# seq_len(length(wavefrm)), col = 'gray', lty = 2)
}
if (type == "ac") {
ac <- acf(x = wavefrm, lag.max = length(wavefrm), type = "covariance",
demean = FALSE, plot = FALSE)
perd <- ac$acf/max(ac$acf)
}
tpks <- seewave::fpeaks(cbind(seq_len(length(perd)), perd), plot = FALSE,
threshold = 0.5)
if (nrow(tpks) > npeak) {
tpks <- tpks[1:npeak, ]
}
segment_df <- data.frame(selec = seq_len(nrow(tpks)), pos = tpks[,
1], peak = tpks[, 2])
# get mean number of sample between peaks
mean_dist_peak <- round(mean(diff(segment_df$pos)))
segment_df$start <- segment_df$pos - mean_dist_peak/2
segment_df$end <- segment_df$pos + mean_dist_peak/2
# fix if values are out of wavefrm size
if (segment_df$start[1] > 0) {
segment_df$start[1] <- 0
}
if (segment_df$end[nrow(segment_df)] > length(wavefrm)) {
segment_df$end[nrow(segment_df)] <- length(wavefrm)
}
# extract segments into a list
segments <- lapply(seq_len(nrow(segment_df)), function(x) {
wavefrm[segment_df$start[x]:segment_df$end[x]]
})
# make all the same number of samples
segments <- lapply(segments, function(x) {
approx(x, n = max(sapply(segments, length)))$y
})
# normalize between 1, -1
segments <- lapply(segments, function(x) {
x/max(x)
})
# put all segments in a data frame
segments <- as.data.frame(segments, col.names = seq_len(length(segments)))
# compute mean segment
mean_segment <- rowMeans(segments)
if (plot) {
mean_segment_df <- data.frame(time = seq(0, 1, length.out = nrow(segments)),
mean.amp = rowMeans(segments), sd.amp = apply(segments,
1, sd))
gg <- ggplot(data = mean_segment_df, mapping = aes(x = time,
y = mean.amp)) + geom_line(color = wave_col) + geom_ribbon(aes(ymin = mean.amp -
sd.amp, ymax = mean.amp + sd.amp), alpha = 0.2) + theme_classic(base_size = 25)
print(gg)
}
if (mean) {
return(mean_segment)
} else {
return(segments)
}
}4.1 Test function
mean period +/- standard deviation using autocorrelation
4.2 Get mean schroeder cycles
4.2.1 EMD method
Code
est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")
mean_schroeders <- warbleR:::pblapply_wrblr_int(cl = 20, seq_len(nrow(est_schr)),
function(x) {
wave <- read_wave(est_schr, index = x)
seg <- try_na(mean_segment(wave, plot = FALSE, mean = FALSE,
type = "EMD"))
return(seg)
})
names(mean_schroeders) <- paste(est_schr$f0, est_schr$label, sep = "-")
saveRDS(mean_schroeders, "./data/processed/mean_random_start_schroeders_emd.RDS")Code
mean_schroeders <- readRDS("./data/processed/mean_random_start_schroeders_emd.RDS")
mean_schroeders <- mean_schroeders[!sapply(mean_schroeders, function(x) is.na(x[[1]][1]))]
mean_schroeders_list <- lapply(seq_len(length(mean_schroeders)), function(x) {
data.frame(schroeder = names(mean_schroeders)[x], time = seq(0,
1, length.out = nrow(mean_schroeders[[x]])), mean.amp = rowMeans(mean_schroeders[[x]]),
sd.amp = apply(mean_schroeders[[x]], 1, sd))
})
mean_schroeders_df <- do.call(rbind, mean_schroeders_list)
ggplot(data = mean_schroeders_df[mean_schroeders_df$schroeder %in%
unique(mean_schroeders_df$schroeder)[1:40], ], mapping = aes(x = time,
y = mean.amp)) + geom_line(color = wave_col) + geom_ribbon(aes(ymin = mean.amp -
sd.amp, ymax = mean.amp + sd.amp), alpha = 0.2) + theme_classic(base_size = 5) +
facet_wrap("~ schroeder", ncol = 5, scales = "free_y")
4.2.2 Autocorrelation method
Code
est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")
mean_schroeders <- warbleR:::pblapply_wrblr_int(cl = 20, seq_len(nrow(est_schr)),
function(x) {
wave <- read_wave(est_schr, index = x)
seg <- try_na(mean_segment(wave, plot = FALSE, mean = FALSE,
type = "ac", thinning = 0.8))
return(seg)
})
names(mean_schroeders) <- paste(est_schr$f0, est_schr$label, sep = "-")
saveRDS(mean_schroeders, "./data/processed/mean_random_start_schroeders_ac.RDS")Code
mean_schroeders <- readRDS("./data/processed/mean_random_start_schroeders_ac.RDS")
mean_schroeders <- mean_schroeders[!sapply(mean_schroeders, function(x) is.na(x[[1]][1]))]
mean_schroeders_list <- lapply(seq_len(length(mean_schroeders)), function(x) {
data.frame(schroeder = names(mean_schroeders)[x], time = seq(0,
1, length.out = nrow(mean_schroeders[[x]])), mean.amp = rowMeans(mean_schroeders[[x]]),
sd.amp = apply(mean_schroeders[[x]], 1, sd))
})
mean_schroeders_df <- do.call(rbind, mean_schroeders_list)
ggplot(data = mean_schroeders_df[mean_schroeders_df$schroeder %in%
unique(mean_schroeders_df$schroeder)[1:40], ], mapping = aes(x = time,
y = mean.amp)) + geom_line(color = wave_col) + geom_ribbon(aes(ymin = mean.amp -
sd.amp, ymax = mean.amp + sd.amp), alpha = 0.2) + theme_classic(base_size = 5) +
facet_wrap("~ schroeder", ncol = 5, scales = "free_y")
5 Measuring schroeder dissimilarity
5.1 Dynamic-time warping pairwise distance
- Both Schroeders have the same length
- One is duplicated and the other one is slide across the duplicated one
- The minimum DTW distance is kept as a dissimilarity measure
Code
mean_schroeders <- readRDS("./data/processed/mean_random_start_schroeders_ac.RDS")
nms <- names(mean_schroeders)
# nms <- grep(pattern = '200|400|600|800', nms, value = TRUE)
cmbs <- t(combn(nms, 2))
min_dist_l <- pbapply::pbsapply(cl = 22, 1:nrow(cmbs), function(x) {
s1 <- rowMeans(mean_schroeders[[cmbs[x, 1]]])
s2 <- rowMeans(mean_schroeders[[cmbs[x, 2]]])
# make same length if (length(s1) != length(s2))
s1 <- approx(s1, n = 100)$y
s2 <- approx(s2, n = 100)$y
# duplicate 1
s1 <- rep(s1, 2)
# run dtw over longer vector
dists <- vapply(seq_len(length(s1) - length(s2)), function(x) {
segment <- s1[x:min(c(x + length(s2) - 1), length(s1))]
dtw_dist <- warbleR::try_na(dtw::dtwDist(mx = rbind(s2, segment)))
return(dtw_dist[1, 2])
}, FUN.VALUE = numeric(1))
return(data.frame(schr1 = cmbs[x, 1], schdr2 = cmbs[x, 2], min(dists)))
})
min_dists <- do.call(rbind, min_dist_l)
min_dists <- as.data.frame(matrix(min_dists[, 1], ncol = 3, byrow = TRUE))
names(min_dists) <- c("schr1", "schr2", "dist")
min_dists$dist <- as.numeric(min_dists$dist)
saveRDS(min_dists, "./data/processed/dtw_distance_ac_segments_random_start.RDS")6 Compare dissimilarity between schroeders using different methods
We try 10 methods for measuring acoustic structure:
- Dynamic time warping (DTW) distances from schroeder waveform
- Raven’s spectrographic features
- warbleR’s spectrographic features
- warbleR’s MFCC features
- warbleR’s cross-correlation of fourier spectrograms
- Raven’s cross-correlation of fourier spectrograms
- Raven’s waveform correlation
- warbleR’s cross-correlation of mel-frequency spectrograms
- DTW on warbleR’s dominant frequency contours
- DTW on Raven’s peak frequency contours
Statistical modeling: - Multiple Regression on distance Matrices - Model:
\[\begin{align*}
Dissimilarity &\sim frequency + components + sign
\end{align*}\] - Response values scaled to make effect sizes comparable across models - Predictors were coded as pairwise binary matrices in which 0 means that calls in a dyad belong to the same level and 1 means calls belong to different levels
6.1 Waveform DTW
Code
min_dists <- readRDS("./data/processed/dtw_distance_ac_segments_random_start.RDS")
min_dists$dist <- scale(min_dists$dist)
dist_tri <- PhenotypeSpace::rectangular_to_triangular(min_dists)
freq_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
"-"), "[[", 1)))
comp_bi_tri <- as.dist(binary_triangular_matrix(group = gsub("_n|_p",
"", sapply(strsplit(rownames(dist_tri), "-"), "[[", 2))))
sign_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
"_"), "[[", 2)))
rect_var <- cbind(as.dist(dist_tri), freq_bi_tri, comp_bi_tri, sign_bi_tri)
colnames(rect_var) <- c("dtw_dist", "frequency", "components", "sign")
dtw_wv_mod <- MRM2(formula = dtw_dist ~ frequency + components + sign,
nperm = 10000, mat = rect_var)
saveRDS(dtw_wv_mod, "./data/processed/matrix_correlation_dtw_distance_random_start.RDS")6.2 Raven’s spectrographic features
Code
rav_dat <- imp_raven(path = "./data/processed", files = "random_start_schroeder_master_raven_measurements.txt",
warbler.format = TRUE, all.data = TRUE)
rav_dat <- rav_dat[grep(pattern = "marker", rav_dat$f0, invert = TRUE),
]
rav_dat <- rav_dat[, grep("contour|PFC Slope", ignore.case = TRUE,
x = names(rav_dat), invert = TRUE)]
pca <- prcomp(x = rav_dat[, names(rav_dat) %in% c("Agg Entropy (bits)",
"Avg Entropy (bits)", "Avg Power Density (dB FS/Hz)", "BW 50% (Hz)",
"BW 90% (Hz)", "Center Freq (Hz)", "Center Time Rel.", "Delta Freq (Hz)",
"Dur 50% (s)", "Dur 90% (s)", "Energy (dB FS)", "Freq 25% (Hz)",
"Freq 5% (Hz)", "Freq 75% (Hz)", "Freq 95% (Hz)", "Inband Power (dB FS)",
"Max Entropy (bits)", "Max Freq (Hz)", "Min Entropy (bits)", "Peak Freq (Hz)",
"PFC Avg Slope (Hz/ms)", "PFC Max Freq (Hz)", "PFC Max Slope (Hz/ms)",
"PFC Min Freq (Hz)", "PFC Min Slope (Hz/ms)", "PFC Num Inf Pts",
"Peak Power Density (dB FS/Hz)", "Peak Time (s)", "Peak Time Relative",
"Time 25% (s)", "Time 25% Rel.", "Time 5% (s)", "Time 5% Rel.",
"Time 75% (s)", "Time 75% Rel.", "Time 95% (s)", "Time 95% Rel.")],
scale. = TRUE)
rav_dat$pc1 <- pca$x[, 1]
rav_dat$pc2 <- pca$x[, 2]
rav_dat$comp <- sapply(strsplit(rav_dat$label, "_"), "[[", 1)
rav_dat$sign <- sapply(strsplit(rav_dat$label, "_"), "[[", 2)
rav_dat$label <- paste(rav_dat$f0, rav_dat$label, sep = "-")
dist_tri <- dist(rav_dat[, c("pc1", "pc2")])
freq_bi_tri <- as.dist(binary_triangular_matrix(group = rav_dat$f0))
comp_bi_tri <- as.dist(binary_triangular_matrix(group = rav_dat$comp))
sign_bi_tri <- as.dist(binary_triangular_matrix(group = rav_dat$sign))
rect_var <- cbind(as.dist(scale(as.matrix(dist_tri))), freq_bi_tri,
comp_bi_tri, sign_bi_tri)
colnames(rect_var) <- c("rav_dist", "frequency", "components", "sign")
rav_mod <- MRM2(formula = rav_dist ~ frequency + components + sign,
nperm = 10000, mat = rect_var)
saveRDS(rav_mod, "./data/processed/matrix_correlation_raven_measurements_distance_random_start.RDS")Code
$coef
rav_dist pval
Int -1.308738783 0.0001
frequency 0.879062058 0.0001
components 0.655162497 0.0001
sign -0.004605695 0.5591
$r.squared
R2 pval
0.09958528 0.00010000
$F.test
F F.pval
1587.7273 0.0001 6.3 warbleR’s spectrographic features
Code
est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")
wsp <- spectro_analysis(est_schr, harmonicity = TRUE, nharmonics = 5,
parallel = 20)
# keep columns with no NAs
wsp <- wsp[complete.cases(wsp), ]
pca <- prcomp(x = wsp[, -c(1:3)])
wsp$pc1 <- pca$x[, 1]
wsp$pc2 <- pca$x[, 2]
wsp_tri <- dist(wsp[, c("pc1", "pc2")])
freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(wsp$sound.files,
"_"), "[[", 2), start = 4, 6)))
comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(wsp$sound.files,
"_"), "[[", 3), start = 7, 8)))
sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
wsp$sound.files), "pos", "neg")))
rect_var <- cbind(as.dist(scale(as.matrix(wsp_tri))), freq_bi_tri,
comp_bi_tri, sign_bi_tri)
colnames(rect_var) <- c("wrbl_sp", "frequency", "components", "sign")
wrbl_mod <- MRM2(formula = wrbl_sp ~ frequency + components + sign,
nperm = 10000, mat = rect_var)
saveRDS(wrbl_mod, "./data/processed/matrix_correlation_warbler_measurements_distance_random_start.RDS")Code
$coef
wrbl_sp pval
Int -1.310541864 0.0001
frequency 1.536291158 0.0001
components 0.104840203 0.0002
sign -0.005210514 0.5868
$r.squared
R2 pval
0.3557964 0.0001000
$F.test
F F.pval
4116.8767 0.0001 6.4 warbleR’s MFCC features
Code
est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")
wsp <- mfcc_stats(est_schr)
# keep columns with no NAs
wsp <- wsp[complete.cases(wsp), ]
pca <- prcomp(x = wsp[, -c(1:3)])
wsp$pc1 <- pca$x[, 1]
wsp$pc2 <- pca$x[, 2]
wsp_tri <- dist(wsp[, c("pc1", "pc2")])
freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(wsp$sound.files,
"_"), "[[", 2), start = 4, 6)))
comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(wsp$sound.files,
"_"), "[[", 3), start = 7, 8)))
sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
wsp$sound.files), "pos", "neg")))
rect_var <- cbind(as.dist(scale(as.matrix(wsp_tri))), freq_bi_tri,
comp_bi_tri, sign_bi_tri)
colnames(rect_var) <- c("wrbl_sp", "frequency", "components", "sign")
mfcc_wrbl_mod <- MRM2(formula = wrbl_sp ~ frequency + components +
sign, nperm = 10000, mat = rect_var)
saveRDS(mfcc_wrbl_mod, "./data/processed/mfcc_warbler_distance_random_start.RDS")6.5 warbleR’s fourier cross-correlation
Code
xc <- readRDS("./data/processed/fourier_cross_correlation_random_start_schroeder.RDS")
freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(xc),
"_"), "[[", 2), start = 4, 6)))
comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(xc),
"_"), "[[", 3), start = 7, 8)))
sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
rownames(xc)), "pos", "neg")))
rect_var <- cbind(as.dist(scale(1 - xc)), freq_bi_tri, comp_bi_tri,
sign_bi_tri)
colnames(rect_var) <- c("fourier_xc", "frequency", "components", "sign")
xc_mod <- MRM2(formula = fourier_xc ~ frequency + components + sign,
nperm = 10000, mat = rect_var)
saveRDS(xc_mod, "./data/processed/matrix_correlation_fourier_cross_correlation_random_start.RDS")Code
$coef
fourier_xc pval
Int -1.78485496 0.0001
frequency 1.98172165 0.0001
components 0.14076790 0.0001
sign -0.02157745 0.0081
$r.squared
R2 pval
0.5070224 0.0001000
$F.test
F F.pval
14764.6553 0.0001 6.6 Raven’s fourier cross-correlation
Code
xc_rav <- Rraven::imp_corr_mat(file = "random_start_schroeders_raven_fourier_cross_correlation.txt",
path = "./data/processed")[[1]]
freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(xc_rav),
"_"), "[[", 2), start = 4, 6)))
comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(xc_rav),
"_"), "[[", 3), start = 7, 8)))
sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
rownames(xc_rav)), "pos", "neg")))
rect_var <- cbind(as.dist(scale(1 - xc_rav)), freq_bi_tri, comp_bi_tri,
sign_bi_tri)
colnames(rect_var) <- c("fourier_xc_rav", "frequency", "components",
"sign")
xc_rav_mod <- MRM2(formula = fourier_xc_rav ~ frequency + components +
sign, nperm = 10000, mat = rect_var)
saveRDS(xc_rav_mod, "./data/processed/matrix_correlation_raven_fourier_cross_correlation_random_start.RDS")Code
$coef
fourier_xc_rav pval
Int -1.7370728 0.0001
frequency 1.7351569 0.0001
components 0.2531801 0.0001
sign 0.0246393 0.0072
$r.squared
R2 pval
0.3757609 0.0001000
$F.test
F F.pval
8700.3890 0.0001 6.7 Raven’s waveform correlation
Code
wav_cr_rav <- Rraven::imp_corr_mat(file = "random_start_schroeders_raven_waveform_correlation.txt",
path = "./data/processed")[[1]]
freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(wav_cr_rav),
"_"), "[[", 2), start = 4, 6)))
comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(wav_cr_rav),
"_"), "[[", 3), start = 7, 8)))
sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
rownames(wav_cr_rav)), "pos", "neg")))
rect_var <- cbind(as.dist(scale(1 - wav_cr_rav)), freq_bi_tri, comp_bi_tri,
sign_bi_tri)
colnames(rect_var) <- c("fourier_wav_cr_rav", "frequency", "components",
"sign")
wav_cr_rav_mod <- MRM2(formula = fourier_wav_cr_rav ~ frequency +
components + sign, nperm = 10000, mat = rect_var)
saveRDS(wav_cr_rav_mod, "./data/processed/matrix_correlation_raven_waveform_correlation_random_start.RDS")Code
$coef
fourier_wav_cr_rav pval
Int -1.70508855 0.0001
frequency 1.88607629 0.0001
components 0.04258164 0.0061
sign 0.11819114 0.0001
$r.squared
R2 pval
0.4953648 0.0001000
$F.test
F F.pval
14188.1430 0.0001 6.8 warbleR’s mel-frequency cross-correlation
Code
xc <- readRDS("./data/processed/mel_frequency_cross_correlation_schroeder_random_start.RDS")
freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(xc),
"_"), "[[", 2), start = 4, 6)))
comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(xc),
"_"), "[[", 3), start = 7, 8)))
sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
rownames(xc)), "pos", "neg")))
rect_var <- cbind(as.dist(scale(1 - xc)), freq_bi_tri, comp_bi_tri,
sign_bi_tri)
colnames(rect_var) <- c("mel_xc", "frequency", "components", "sign")
mxc_mod <- MRM2(formula = mel_xc ~ frequency + components + sign,
nperm = 10000, mat = rect_var)
saveRDS(mxc_mod, "./data/processed/matrix_correlation_mel_cross_correlation_random_start.RDS")Code
$coef
mel_xc pval
Int -1.2613675324 0.0001
frequency 1.2552094577 0.0001
components 0.2486046157 0.0001
sign 0.0008366943 0.9119
$r.squared
R2 pval
0.1794016 0.0001000
$F.test
F F.pval
3138.4764 0.0001 6.9 DTW on warbleR’s dominant frequency contours
Code
source("~/Dropbox/R_package_testing/warbleR/R/freq_DTW.R")
source("~/Dropbox/R_package_testing/warbleR/R/internal_functions.R")
est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")
dtw_dists <- freq_DTW(est_schr, type = "dominant", img = FALSE)
saveRDS(dtw_dists, "./data/processed/dtw_distance_dominant_frequency_schroeder_random_start.RDS")Code
dtw_dists <- readRDS("./data/processed/dtw_distance_dominant_frequency_schroeder_random_start.RDS")
freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(dtw_dists),
"_"), "[[", 2), start = 4, 6)))
comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(dtw_dists),
"_"), "[[", 3), start = 7, 8)))
sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
rownames(dtw_dists)), "pos", "neg")))
rect_var <- cbind(as.dist(scale(dtw_dists)), freq_bi_tri, comp_bi_tri,
sign_bi_tri)
colnames(rect_var) <- c("dtw_dists", "frequency", "components", "sign")
dtw_dom_mod <- MRM2(formula = dtw_dists ~ frequency + components +
sign, nperm = 10000, mat = rect_var)
saveRDS(dtw_dom_mod, "./data/processed/matrix_correlation_dtw_distances_dominant_frequency_contours_random_start.RDS")Code
$coef
dtw_dists pval
Int -0.456171403 0.0001
frequency 0.526558411 0.0001
components 0.157208372 0.0001
sign -0.004413661 0.5422
$r.squared
R2 pval
0.02564736 0.00010000
$F.test
F F.pval
377.8765 0.0001 6.10 DTW on Raven’s peak frequency contours
Code
source("~/Dropbox/R_package_testing/warbleR/R/freq_DTW.R")
source("~/Dropbox/R_package_testing/warbleR/R/internal_functions.R")
rav_dat <- imp_raven(path = "./data/processed", files = "random_start_schroeder_master_raven_measurements.txt",
warbler.format = TRUE, all.data = TRUE)
rav_dat <- rav_dat[grep(pattern = "marker", rav_dat$f0, invert = TRUE),
]
# rav_dat <- rav_dat[, grep('contour|PFC Slope', ignore.case =
# TRUE, x = names(rav_dat), invert = TRUE)]
rav_freq_cntr <- extract_ts(rav_dat, ts.column = "Peak Freq Contour (Hz)",
equal.length = TRUE)
est_schr <- readRDS("./data/processed/extended_selection_table_random_start_schroeders.RDS")
rav_dtw_dists <- freq_DTW(est_schr, ts.df = rav_freq_cntr, type = "dominant",
img = FALSE)
rownames(rav_dtw_dists) <- rav_dat$orig.sound.file
saveRDS(rav_dtw_dists, "./data/processed/dtw_distance_raven_dominant_frequency_random_start_schroeder.RDS")Code
rav_dtw_dists <- readRDS("./data/processed/dtw_distance_raven_dominant_frequency_random_start_schroeder.RDS")
freq_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(rav_dtw_dists),
"_"), "[[", 2), start = 4, 6)))
comp_bi_tri <- as.dist(binary_triangular_matrix(group = substr(sapply(strsplit(rownames(rav_dtw_dists),
"_"), "[[", 3), start = 7, 8)))
sign_bi_tri <- as.dist(binary_triangular_matrix(group = ifelse(grepl("pos",
rownames(rav_dtw_dists)), "pos", "neg")))
rect_var <- cbind(as.dist(scale(rav_dtw_dists)), freq_bi_tri, comp_bi_tri,
sign_bi_tri)
colnames(rect_var) <- c("rav_dtw_dists", "frequency", "components",
"sign")
source("~/Dropbox/Projects/geographic_call_variation_yellow-naped_amazon/scripts/MRM2.R")
dtw_rav_peak_mod <- MRM2(formula = rav_dtw_dists ~ frequency + components +
sign, nperm = 10000, mat = rect_var)
saveRDS(dtw_rav_peak_mod, "./data/processed/matrix_correlation_dtw_distances_raven_dominant_frequency_contours_random_start.RDS")Code
$coef
rav_dtw_dists pval
Int -0.8465155615 0.0001
frequency 1.1803304963 0.0001
components 0.0468889213 0.0323
sign 0.0009260393 0.9175
$r.squared
R2 pval
0.1363664 0.0001000
$F.test
F F.pval
2266.7382 0.0001 6.11 Combined results
Code
mods <- list(dtw_wv_mod = dtw_wv_mod, wrbl_mod = wrbl_mod, mfcc_wrbl_mod = mfcc_wrbl_mod,
rav_mod = rav_mod, xc_mod = xc_mod, mxc_mod = mxc_mod, dtw_dom_mod = dtw_dom_mod,
dtw_rav_peak_mod = dtw_rav_peak_mod, wav_cr_rav_mod = wav_cr_rav_mod,
xc_rav_mod = xc_rav_mod)
names(mods) <- c("DTW waveform period", "warbleR's spectro features",
"warbleR's MFCC features", "Raven's spectro features", "warbleR's Fourier cross-correlation",
"Mel cross-correlation", "DTW warbleR's dominant freq", "DTW Raven's peak freq",
"waveform correlation", "Raven's Fourier cross-correlation")
estimates <- do.call(rbind, lapply(seq_along(mods), function(x) {
Y <- data.frame(mods[[x]]$coef[-1, ])
Y$rel_coef <- Y[, 1]/max(Y[, 1])
Y$mod <- names(mods)[x]
Y$predictor <- rownames(Y)
names(Y) <- c("coef", "p", "rel_coef", "model", "predictor")
return(Y)
}))
estimates$rel_coef <- ifelse(estimates$p < 0.05, estimates$rel_coef,
0)
estimates$signif <- ifelse(estimates$p < 0.05, "p < 0.05", "p >= 0.05")
ggplot(estimates, aes(x = predictor, y = model, fill = rel_coef)) +
geom_tile() + coord_equal() + scale_fill_gradient2(low = viridis(10)[3],
high = viridis(10)[7], guide = "none") + geom_text(aes(label = round(coef,
3), color = signif)) + scale_color_manual(values = c("black",
"gray")) + labs(x = "", y = "", color = "P value") + theme_classic() +
theme(axis.text.x = element_text(color = "black", size = 11, angle = 30,
vjust = 0.8, hjust = 0.8))
7 Takeaways
- Amplitude autocorrelation works better at getting periodicity
- Dynamic time warping dissimilarity on mean periods, Raven’s waveform correlation and (to some extent) Fourier cross-correlation did capture sign differences in Schroeders
- Dynamic time warping dissimilarity on mean periods remove differences in frequency when forcing all periods to have the same number of values (i.e. same length)
Session information
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
locale:
[1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
[3] LC_TIME=es_CR.UTF-8 LC_COLLATE=pt_BR.UTF-8
[5] LC_MONETARY=es_CR.UTF-8 LC_MESSAGES=pt_BR.UTF-8
[7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] animation_2.6 numform_0.7.0 ecodist_2.0.9
[4] PhenotypeSpace_0.1.0 ggplot2_3.4.2 baRulho_1.1.0
[7] ohun_1.0.0 Rraven_1.0.13 viridis_0.6.3
[10] viridisLite_0.4.2 warbleR_1.1.28 NatureSounds_1.0.4
[13] tuneR_1.4.4 seewave_2.2.1 formatR_1.11
[16] knitr_1.43 kableExtra_1.3.4
loaded via a namespace (and not attached):
[1] colorspace_2.1-0 rjson_0.2.21 deldir_0.2-10
[4] class_7.3-22 xaringanExtra_0.7.0 rstudioapi_0.14
[7] proxy_0.4-27 spatstat.data_2.1-0 farver_2.1.1
[10] Deriv_4.1.3 remotes_2.4.2 fansi_1.0.4
[13] xml2_1.3.5 codetools_0.2-18 splines_4.1.0
[16] polyclip_1.10-0 jsonlite_1.8.7 vdiffr_1.0.5
[19] packrat_0.9.0 cluster_2.1.2 png_0.1-7
[22] spatstat.sparse_2.0-0 compiler_4.1.0 httr_1.4.4
[25] backports_1.4.1 assertthat_0.2.1 Matrix_1.5-4.1
[28] fastmap_1.1.1 cli_3.6.1 htmltools_0.5.5
[31] tools_4.1.0 igraph_1.5.0.1 gtable_0.3.3
[34] glue_1.6.2 dplyr_1.0.10 Rcpp_1.0.11
[37] raster_3.4-13 vctrs_0.6.3 svglite_2.1.0
[40] nlme_3.1-162 xfun_0.39 stringr_1.5.0
[43] brio_1.1.3 testthat_3.1.10 rvest_1.0.3
[46] lifecycle_1.0.3 goftest_1.2-2 MASS_7.3-60
[49] scales_1.2.1 spatstat.core_2.3-0 spatstat.utils_2.2-0
[52] parallel_4.1.0 Sim.DiffProc_4.8 yaml_2.3.7
[55] pbapply_1.7-2 gridExtra_2.3 rpart_4.1-15
[58] stringi_1.7.12 e1071_1.7-13 checkmate_2.2.0
[61] permute_0.9-5 rlang_1.1.1 pkgconfig_2.0.3
[64] systemfonts_1.0.4 dtw_1.23-1 bitops_1.0-7
[67] evaluate_0.21 lattice_0.21-8 tensor_1.5
[70] sf_1.0-14 labeling_0.4.2 htmlwidgets_1.5.4
[73] tidyselect_1.2.0 magrittr_2.0.3 R6_2.5.1
[76] fftw_1.0-7 generics_0.1.3 DBI_1.1.3
[79] pillar_1.9.0 withr_2.5.0 mgcv_1.8-42
[82] units_0.8-2 abind_1.4-5 RCurl_1.98-1.12
[85] sp_1.5-1 tibble_3.2.1 crayon_1.5.2
[88] KernSmooth_2.23-21 utf8_1.2.3 spatstat.geom_2.2-2
[91] rmarkdown_2.23 grid_4.1.0 sketchy_1.0.2
[94] vegan_2.5-7 digest_0.6.33 classInt_0.4-9
[97] webshot_0.5.4 signal_0.7-7 munsell_0.5.0